%% Gaussian Processes
% Tobias Siegfried, 06/03/2008

curP = pwd;
cd('/Users/tobiassiegfried/Documents/Science/matlab/agents/RLGP');
clear all
close all
%% 1D
figure(1)
m = [2, 3, 5, 10, 20];
for i = 1:5
    n = m(i);
    rand('state',18);
    randn('state',20);
    covfunc = {'covSum', {'covSEiso','covNoise'}};
    loghyper = [log(1.0); log(1.0); log(0.1)];
    x = 15*(rand(n,1)-0.5);
    y = chol(feval(covfunc{:}, loghyper, x))'*randn(n,1);
    plot(x, y, 'k+', 'MarkerSize', 17)
    %loghyper = minimize([-1; -1; -1], 'gpr', -100, covfunc, x, y)
    %keyboard
    % compute
    xstar = linspace(-7.5,7.5,201)'; % 201 test points / discretization
    [mu S2] = gpr(loghyper, covfunc, x, y, xstar) % predictions
    S2 = S2 - exp(2*loghyper(3));
    clf
    f = [mu+2*sqrt(S2);flipdim(mu-2*sqrt(S2),1)]; % standard noise-free pointwise errorbars
    fill([xstar; flipdim(xstar,1)], f, [7 7 7]/8, 'EdgeColor', [7 7 7]/8);
    hold on
    plot(xstar,mu,'b-','LineWidth',1);
    plot(x, y, 'k+', 'MarkerSize', 17);
    axis([-8 8 -3 3]);
    keyboard
end
%% 2D
load data_6darm.mat
mean(X), std(X), mean(y), std(y)
offset = mean(y);
y = y - offset;         % centre targets around 0.
covfunc = {'covSum', {'covSEard','covNoise'}};
logtheta0 = [0; 0; 0; 0; 0; 0; 0; log(sqrt(0.1))];
[logtheta, fvals, iter] = minimize(logtheta0, 'gpr', -10, covfunc, X, y);
exp(logtheta)
[fstar S2] = gpr(logtheta, covfunc, X, y, Xstar);
fstar = fstar + offset;  % add back offset to get true prediction
res = fstar-ystar;  % residuals
mse = mean(res.^2)
pll = -0.5*mean(log(2*pi*S2)+res.^2./S2)
subplot(211), plot(res,'.'), ylabel('residuals'), xlabel('test case')
subplot(212), plot(sqrt(S2),'.'), ylabel('predictive std deviation'), xlabel('test case')
%% go back
cd curP
%% test - 
figure(11)
xes = [-2:.2:2];
[xo,y] = meshgrid(xes);
xstarT = (-2:.2:2)'; xstar1 = repmat(xstarT,length(xstarT),1);
xstar2 = xstarT'; xstar2 = repmat(xstar2,length(xstarT),1);
xstar2 = reshape(xstar2,size(xstar2,1)*size(xstar2,2),1);
xstar = [xstar1 xstar2];
z = xo.*exp(-xo.^2-y.^2); %*randn(length(x),length(y));
surfc(xo,y,z)
% pick some samples randomly
% threshold
t = 0.5;
mRand = rand(length(xo),length(y));
mRand(mRand>t) = 1;
mRand(mRand~=1) = 0;
surfc(xo,y,z.*mRand)
[x1,x2] = find(mRand);
x = [xes(x1)' xes(x2)'];
xi = sub2ind(size(mRand),x1,x2);
yRes = z(xi);
%%
figure(1)
n = sum(mRand(:));
rand('state',18);
randn('state',20);
covfunc = {'covSum', {'covSEiso','covNoise'}};
loghyper = [log(1.0); log(1.0); log(0.1)];
%x = 15*(rand(n,1)-0.5);
%y = chol(feval(covfunc{:}, loghyper, x))'*randn(n,1);
stem3(x(:,1), x(:,2), yRes,'+')
%loghyper = minimize([-1; -1; -1], 'gpr', -100, covfunc, x, y)
%keyboard
%% compute
%xstar = linspace(-7.5,7.5,201)'; % 201 test points / discretization
figure(2)
[mu S2] = gpr(loghyper, covfunc, x, yRes, xstar); % predictions
S2 = S2;% - exp(2*loghyper(3));
clf
f = [mu+2*sqrt(S2);flipdim(mu-2*sqrt(S2),1)]; % standard noise-free pointwise errorbars
fill([xstar; flipdim(xstar,1)], f, [7 7 7]/8, 'EdgeColor', [7 7 7]/8);
hold on
plot(xstar,mu,'b-','LineWidth',1);
plot(x, yRes, 'k+', 'MarkerSize', 17);
axis([-8 8 -3 3]);
%%
figure(4)
surfc(xo,y,reshape(mu,length(xo),length(xo))')
hold on
stem3(x(:,1), x(:,2), yRes,'k+')
hold off



